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Matteo Beccaria^, Carlo Presilla*^, Gian Fabrizio De Angelis'^ and Giovanni Jona-Lasinio'^ 

** Dipartiniento di Fisica, Universita di Lecce, Via Arnesano, 73100 Lecce, Italy 

^ Dipartimento di Fisica, Universita di Roma "La Sapienza", Piazzale A. Moro 2, 00185 Roma, Italy 
and INFM Unita di Ricerca di Roma "La Sapienza" 

"^ Dipartimento di Fisica, Universita di Lecce, Via Arnesano, 73100 Lecce, Italy 

'^ Dipartimento di Fisica, Universita di Roma "La Sapienza", Piazzale A. Moro 2, 00185 Roma, Italy 

We describe an exact Feynman-Kac type formula to represent the dynamics of fermionic lattice systems. In 
this approach the real time or Euclidean time dynamics is expressed in terms of the stochastic evolution of a 
collection of Poisson processes. From this formula we derive a family of algorithms for Monte Carlo simulations, 
parametrized by the jump rates of the Poisson processes. 



Quantum Monte Carlo methods are powerful 
techniques for the numerical evaluation of the 
properties of quantum lattice systems. In the case 
of fermion systems |l|-Q] there are special features 
connected with the anticommutativity of the vari- 
ables involved. In a recent paper H progress 
has been made by providing an exact probabilis- 
tic representation for the dynamics of a Hubbard 
model. Here we illustrate the basic formula while 
for details we refer to H. 

Let us consider the Hubbard Hanriltonian 

|A| |A| 

i=l j=i+l (T=tX 
|A| 

+ IZt' ^^T^n c.Va' (1) 

where A C Z"^ is a finite d-dimensional lattice 
with cardinality |A|, {1,...,|A|} some total or- 
dering of the lattice points, and Cia the usual an- 
ticommuting destruction operators at site i and 
spin index a. In this paper, we are interested in 
evaluating the matrix elements (n'|e~^*|n) where 
n = (ni|, ni|, . . . , ni\n, '^lAli) ^^^ the occupation 
numbers taking the values or 1. The total num- 
ber of fermions per spin component is a conserved 
quantity, therefore we consider only configura- 



tor a =ti- In the following we shall use the mod 
2 addition n ® n' = {n + n') mod 2. 

Let P = {(i,j),l < I < j < \A\ : rfc ^ 0} 
and |r| its cardinality. For simplicity, we will as- 
sume that rjij = ?7 if (i,j) € P and ji = 7. By 
introducing 
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where 1^^ = (0, . . . , 0, Ua, 0, . . . , 0), and 
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ni]n.,i, 



the following representation holds 

(n'|e-^*|n) = E (j„,^„eN'A^*) 

7W* = expi Y^ J2 h^S ["^P^^ 
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X Ay^(n©N^)]diVfj„ 

- / V{n®W)ds + 2\T\pt\. (5) 

Here, {Nf^J, (ij) £ P, is a family of 2|P| in- 
dependent Poisson processes with paranreter p 



and N* = (iV*^,7Vf^,...,7V*^|^,iV*^|^) are 2|A| 
stochastic processes defined as 



NL = 



Nl 
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We remind that a Poisson process iV* with pa- 
rameter p is a jump process characterized by the 
foUowing probabihties: 
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(7) 



Its trajectories are piecewise-constant increasing 
integer- valued functions continuous from the left. 
The stochastic integral J dN* is just an ordinary 
Stieltjes integral 



lo.t) 



f{s,N')dN' 



:: Sk<t 



where Sk are random jump times having proba- 
bility density p{s) = pe~P'^. Finally, the symbol 
E(. . .) is the expectation of the stochastic func- 
tional within braces. We emphasize that a similar 
representation holds for the real time matrix ele- 
ments (n'|e~'"^*|n). 

Summarizing, we associate to each rjij 7^ a 
link connecting the sites i and j and assign to it a 
pair of Poisson processes N-j^ with a =TJ.- Then, 
we assign to each site i and spin component a 
a stochastic process iV|^ which is the sum of all 
the processes associated with the links incoming 
at that site and having the same spin component. 
A jump in the link process Nf,^ implies a jump in 
both the site processes N^^ and Nj^. Equations 
(^) and (||) are immediately generalizable to non 
identical parameters rjij and %. In this case, it 
may be convenient to use Poisson processes iV*^^ 
with different parameters pij^- 

In order to construct an efficient algorithm for 
evaluating (Q^), we start by observing that the 
functions Aycr(n © N*) vanish when the occupa- 
tion numbers riia ® Nf^ and rija © N^^ are equal. 
We say that for a given value of a the link ij is 
active at time s if Aijcr(n® N") 7^ 0. We shall see 
in a moment that only active links are relevant. 
Let us consider how the stochastic integral in (^) 
builds up along a trajectory defined by consider- 
ing the time ordered succession of jumps in the 



family {N^,^}. The contribution to the stochastic 
integral in the exponent of (|5|) at the first jump 
time of a link, for definiteness suppose that the 
link iiji with spin component ci jumps first at 
time si, is 

log [r]p-^\^,j,^, (n e N'*! )] 9{t - si), 

where N'*^ — due the assumed left continuity. 
Therefore, if the link iijiai was active at time 
we obtain a finite contribution to the stochas- 
tic integral otherwise we obtain —00. If si > i 
we have no contribution to the stochastic inte- 
gral from this trajectory. If si < t a second jump 
of a link, suppose 12 j2 with spin component CT2, 
can take place at time S2 > si and we obtain a 
contribution 
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(n®N^^)] e{t-S2). 



The analysis can be repeated by considering an 
arbitrary number of jumps. Of course, when the 
stochastic integral is —00, which is the case when 
some A = 0, there is no contribution to the expec- 
tation. The other integral in (^ is an ordinary 
integral of a piecewise constant bounded function. 
We now describe the algorithm. The only tra- 
jectories to be considered are those associated to 
jumps of the active links. This guarantees con- 
servation of the total number of fermions per 
spin component. The active links can be de- 
termined after each jump by inspecting the oc- 
cupation numbers of the sites in the set F ac- 
cording to the rule that the link ij is active for 
the spin component a if riia- + rija ~ 1. We 
start by determining the active links in the ini- 
tial configuration n assigned at time and make 
an extraction with uniform distribution to de- 
cide which of them jumps first. We then extract 
the jump time si according to the probability 
density PAi{s) = Aipexp(— Aips) where Ax is 
the number of active links before the first jump 
takes place. The contribution to yVf* at the time 
of the first jump is therefore, up to the factor 
exp(-2|F|pi), 
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where exp[— (2|r| — Ai)ps] is the probabihty that 
the 2|r| — Ai non active hnks do not jump in 
the time interval s. The contribution of a given 
trajectory is obtained by multiplying the factors 
corresponding to the different jumps until the last 
jump takes place later than t. For a given trajec- 
tory we thus have 

k>l 

^^[Akp-V{n(B'N''k)](s^-Sk_i) q(, _ ^\ 

_j_g[A,p-V(„®N')](t-.,_i) 0(s^, _t)l. (8) 

Here, Ak = A{n © 'N"'') is the number of ac- 
tive links in the interval {sk-i,Sk] and sq = 

0. Note that the exponentially increasing factor 
exp (2|r|pi) in (^) cancels out in the final expres- 
sion of A4*. The analogous expression of A4* for 
real times is simply obtained by replacing rj —> irj 
and 7 — > 17. The algorithm can be improved by 
the usual methods of reconfiguration and impor- 
tance sampling. 

In principle, the algorithms parametrized by p 
are all equivalent as (|^||) holds for any choice 
of the Poisson rates. However, since we estimate 
the expectation values with a finite number of 
trajectories, this may introduce a systematic er- 
ror. It can be shown that the best performance 
is obtained for the natural choice p ^ -q indepen- 
dently of the interaction strength. In this case 
our algorithm coincides with the Green function 
Monte Carlo method in the limit when the latter 
becomes exact as discussed in [||. 
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